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Abstract 

We develop a variational method of deriving stochastic partial differential equations whose solutions follow 
the flow of a stochastic vector field. As an example in one spatial dimension we numerically simulate singular 
solutions (peakons) of the stochastically perturbed Camassa-Holm (CH) equation derived using this method. 
These numerical simulations show that peakon soliton solutions of the stochastically perturbed CH equation 
persist and provide an interesting laboratory for investigating the sensitivity and accuracy of adding stochasticity 
to finite dimensional solutions of stochastic partial differential equations (SPDE). In particular, some choices of 
stochastic perturbations of the peakon dynamics by Wiener noise (canonical Hamiltonian stochastic deformations, 
or CH-SD) allow peakons to interpenetrate and exchange order on the real line in overtaking collisions, although 
this behaviour does not occur for other choices of stochastic perturbations which preserve the Euler-Poincare 
structure of the CH equation (parametric stochastic deformations, or P-SD), and it also does not occur for 
peakon solutions of the unperturbed deterministic CH equation. The discussion raises issues about the science 
of stochastic deformations of finite-dimensional approximations of evolutionary PDE and the sensitivity of the 
resulting solutions to the choices made in stochastic modelling. 
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1 Introduction 

Two main approaches have arisen recently for implementing variational principles in stochastic geometric mechanics. 
In one of them, sometimes called stochastic deformation [2] , the Lagrangian in Hamilton’s principle is the classical 
one, but it is evaluated on underlying stochastic processes and their mean derivatives. This perspective was initially 
motivated by the quantization of classical systems I3IM1I37] and a probabilistic version of Feynman’s path integral 
approach to quantum mechanics. The stochastic deformation approach brings to bear the full theory of stochastic 
partial differential equations (SPDE). For more details about the history of the applications of this approach, in 
particular for fluid dynamics, see [aim- 

Here we will advocate a simpler, and more restricted approach. Our approach is based on a generalisation in m 
of earlier work by Bismut |1] , Lazaro-Cami and Ortega [25] , and Bou-Rabee and Owhadi [5] for stochastic ordinary 
differential equations (SDE), which unifies their Hamiltonian and Lagrangian approaches to temporal stochastic 
dynamics, and extends them to stochastic partial differential equations (SPDE) in the case of cylindrical noise 
in which the spatial dependence is parametric, while temporal dependence is stochastic. The advantage of this 
approach is that the parametric spatial dependence of the stochastic dynamical variables in the resulting SPDEs 
allows essentially finite-dimensional stochastic methods to be applied at each point of space. This feature allows us 
to safely assume from the onset that all the objects we introduce in this context are semimartingales. To distinguish 
between the two approaches, we will call the approach taken here parametric stochastic deformation (P-SD). 

Parametric stochastic partial differential equations (P-SPDE) result in the present approach by applying P-SD to 
a deterministic variational principle. These P-SPDE contain a type of multiplicative, cylindrical, Stratonovich noise 
that depends on both the solution variables and their spatial gradients. This unfamiliar feature does not interfere 
with the passage to the Ito representation, though, since the space variable is treated merely as a parameter when 
dealing with cylindrical noise. That is, one may regard the cylindrical noise process as a finite dimensional stochastic 
process parametrized by x (the space variable). Then, the Stratonovich equation makes analytical sense pointwise, 
for each fixed x. Once this is agreed, then the transformation to ltd by the standard method also makes sense 
pointwise in space. For more details about P-SPDE explained in a fluid dynamics context, see m- 

In this paper we develop an approach for inserting parametric stochastic deformation with cylindrical noise 
into systems of evolutionary partial differential equations which derive from deterministic variational principles 
that are invariant under a Lie group action. The corresponding deterministic dynamical systems are called Euler- 
Poincare equations. The set of Euler-Poincare equations includes the equations of ideal fluid dynamics, which follow 
from variational principles whose Lagrangians satisfy certain invariance properties under smooth invertible maps 
(diffeomorphisms) [271IT5] . 


Objectives. This paper has two main objectives. The first objective is the inclusion of parametric stochastic 
deformation (P-SD) in the variational principle for the EPDiff partial differential equation^ The second objective 
is the numerical study of the statistical effects of parametric and canonically Hamiltonian stochastic deformations 
(CH-SD) on the soliton-like solutions of EPDiff which arise in the deterministic case in one spatial dimension, when 
the Lagrangian in Hamilton’s principle is a Sobolev norm on the continuous vector fields. When the norm is 
chosen this results in the CH equation [8] . We study the stochastic generalised Camassa-Holm equation because its 
singular momentum map m persists under the stochastic deformations we introduce here and thereby allows its 
solutions to be investigated on a finite dimensional invariant manifold. 

^EPDiff is the PDE which arises when the Lagrangian in Hamilton’s principle is a functional of continuous Eulerian vector fields, 
whose flows are smooth invertible maps (diffeomorphisms). EPDiff is the Euler-Poincare equation arising for Lagrangians which are 
invariant under the diffeomorphism group. When the Lagrangian is chosen to be the norm of the vector fields, the EPDiff equation 
becomes the Camassa-Holm equation [^, which is a completely integrable Hamiltonian system and possesses singular soliton solutions 
known for their peaked shape as peakons. 
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1.1 Stochastic EPDiff variational principle 


The action integral for the stochastic variational principle we shall study for EPDiff is a stochastically constrained 
variational principle SS = 0, with action integral, S, given by 


S(u,p,q)= J (^£{u)dt+(p, dq +£uq 


( 1 . 1 ) 


where £{u) is the unperturbed deterministic Lagrangian, written as a functional of velocity vector field u G 3£(K^). 
The angle brackets 

{p,q)v:=J < p{x,t),q{x,t) > dx ( 1 . 2 ) 

denote the spatial integral over the domain of flow of the pairing < p, q > between stochastic dynamical variables 


q, which take values in a tensor space V, and their dual elements p taking values inV*. In (1.1), the quantity p gV 


is a Lagrange multiplier and £uq is the stochastic Lie differential of the dynamical variable q G V with respect to 
a stochastic vector field u(x,t) which is defined by the following sum of a drift velocity u{x,t) and Stratonovich 
stochastic process with cylindrical noise parameterised by spatial position x, 


u{x, t) = u{x, t) dt — ^ ^i{x) o dWi{t). 


(1.3) 


We give a precise definition of the stochastic Lie differential £uq in formula (1.7) in Section 1.2 One may inter¬ 


pret equation (1.3) as the decomposition of a vector field u{x,t) defined at spatial position x and time t into a 


time-dependent drift velocity u{x, f) and a stochastic vector field with cylindrical noise (note that for notational 
convenience we define m as a stochastic differential with respect to the time variable; we will use the tilde to denote 
objects defined as stochastic differentials throughout the rest of this work). The time-independent vector fields £,i{x) 
with i = 1,2,..., M in the cylindrical stochastic process are usually interpreted as “diffusivities” of the stochastic 
vector field, and the choice of these M quantities must somehow be specified from the physics of the problem to be 
considered. In the present considerations, a natural choice will arise from the singular momentum map admitted 
by the deterministic EPDiff equation HZ]. 

The pairing (• j ')v in the Stochastic Variational Principle (SVP) SS = 0 for action functional S in 
with Lagrange multiplier p G T*V enforces the advection condition that the quantity q G V is preserved along 


the Stratonovich stochastic integral curves of the vector field (1.3) for any tensor space V. Namely, q satisfies the 
advection condition 

dq + £uq = 0 . (1.4) 


The advection condition ( |1.4[ ) for the quantities q G V may be regarded as a stochastic constraint imposed on the 
variational principle (O) via the Lagrange multiplier p. 

At this point, we have introduced parametric Stratonovich stochasticity into the variational principle iHl) 
with dynamical variables (u,q,p) through the constraint that the advected quantities q G V should evolve by 


following the Stratonovich stochastic vector field u{x,t) in equation (1.3). This advection law is formulated as a Lie 
differential with respect to the Stratonovich stochastic vector field acting on a tensor space. We use the Stratonovich 
formulation, so the normal rules of calculus apply. For mathematical discussions of Lie derivatives with respect to 
stochastic vector fields, see, e.g., Eiiiia. 


Plan of the paper. After setting out the general theory of Euler-Poincare evolutionary P-SPDE in (x, t) G x K 
in the remainder of this introductory section, we will specialise to one spatial dimension for (x, t) G K x K and study 
the corresponding soliton-like solution behaviour for the case when the Lagrangian £{u) in the variational principle 
is chosen as a Sobolev norm of the vector field u. 

The objective of the remainder of the paper is to use the Stratonovich stochastic EPDiff Theorem proved below 
to study the effects of introducing this type of stochasticity on the interactions of the peakon solutions of the CH 
equation with parametric stochastic deformation (P-SD). The P-SD of the equations of motion for the singular 
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solutions of stochastic EPDiff in one spatial dimension will be introduced via Hamilton’s principle in Section 
and their effects on the numerical solutions will be studied thereafter in comparison with a more general canonical 
Hamiltonian stochastic deformation (CH-SD) in the sense of [HHS], which includes P-SD but can be more general. 

The Fokker-Planck equations for the probability density evolution associated with P-SD and CH-SD of the 
EPDiff equation will be discussed in Section In Section we discuss the numerical algorithm we use, and in 
Sectionj^we present the results of our numerical studies, including sample paths and mean solutions, the probability 
distribution for crossing of singular solutions on the real line, statistics of the first crossing time, effects of noise 
screening, comparison with other types of noise such as additive noise in the canonical momentum equation and 
the results of convergence tests for our stochastic variational integrator. Section]^ is devoted to a brief summary of 
results and discussion of some open problems for possible future work. 


1.2 Stratonovich stochastic EPDiff equation 

If the drift velocity vector field u{x,t) and the diffusivity vector fields satisfy some standard measurability 


and regularity conditions, then the stochastic vector field (1.3) possesses a pathwise unique stochastic flow Ft^s{x). 
By definition, this flow almost surely satisfies Fg fix) = x and the integral equation 


Ft 2 ,six) - Ft^^six) = f u{Ft^s{x),t)dt-'^ f ^i{Ft^six)) o dWfit), 

Jti , Jti 


(1.5) 


or shortly, in the differential form. 


dFt,s{x) = u{Ft^s{x),t) dt -'^^i{Ft,s{x)) o dWfit). (1.6) 

i 

It can be proved that for fixed t, s this flow is mean-square differentiable with respect to the x argument, and also 
almost surely is a diffeomorphism (see m, m, m, m)- These properties allow us to generalise the differential- 
geometric notion of the Lie derivative, which we do in the following definition and theorem. 

Definition 1 (Stochastic Lie differential). Let q he a smooth tensor field. The stochastic Lie differential T^q is the 
almost surely unique stochastic differential satisfying 

d{F,%q) -.= FfJ^q. (1.7) 

Theorem 2. The stochastic Lie differential £^q is almost surely unique and given by 

£uq = £uqdt-^£^^qo dWfit), (1.8) 


where £„ and £c. are the standard Lie derivatives. 


Proof. The proof is a straightforward generalisation of the standard differential-geometric construction of the Lie 
derivative of tensor fields (see [T]). Whenever necessary, we replace time-differentiation with stochastic differentials 


and use the weak property (1.6) of the flow. We first prove (1.8) when q is a smooth, real-valued function. Then 


{Ff^q){x) = qo Ft^six). Formula (1.8) is proved by calculating the stochastic differential d^Ff^q) using the rules 
of Stratonovich calculus. Next we consider the case when g is a smooth vector field. Let G\ be the smooth flow 
of q. Then for fixed t, s the flow Hx of the vector field Ff^q satisfies Gx o Ft ^ = Ff g o H\- From the mean-square 
differentiability of Ft^s we have mean-square differentiability of both sides with respect to A. Differentiating both 
sides with respect to A, evaluating at A = 0, calculating the stochastic differential with respect to t and comparing 


terms, we obtain formula (1.8) for vector fields. For the case when g is a differential one-form we use the property 

' is the dual pairing between one-forms 


Ff^{q, v) = {Ff^q, Ff^v), where v is an arbitrary smooth vector field and 


and vector fields. Calculating the stochastic differential of both sides and using our already established results for 
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functions and vector fields, we prove (1.8) for differential one-forms. It is now straightforward to complete the proof 
for a general tensor field. Almost sure uniqueness of £^<1 follows from our construction and pathwise uniqueness of 
the flow Ft^s- n 

In order to choose the form of the spatial correlations, or diffusivities, ^i(x) of the cylindrical Stratonovich 
stochasticity in (1.3), we notice that the action integral for the variational principle in (1.1) may be rearranged into 
the equivalent form 


S{u,p,q) = 


u)dt (p ^ dq + £uq 


£{u)dt + (p , dq)y - {poq,u) 


X ) ’ 


(1.9) 


where we define the diamond operation (o) in the expression po q £ X* via the real-valued nondegenerate pairing 
(•,-)^:A*xX—between a vector field rj G X and its dual under pairing p G X as 

{poq, rj)y^ := {p, -£rfq)y . ( 1 . 10 ) 

The diamond operation (o) will be instrumental in deriving the Stratonovich form of the stochastic EPDiff equation 
from the stochastic variational principle for the action integral in (1.1), as stated in the following theorem. 

Theorem 3 (Stratonovich Stochastic EPDiff equation). The parametric Stratonovich stochastic deformation in 
the action S{u,p,q) for the stochastic variational principle 5S = 0 for EPDiff given by 

S{u,p,q)= J (^£{u) - {po q , u{x,t))^^dt+ j {p{x,t), dq{x,t))y + j '^{p o q , ^i{x)) dW,{t), 


Lebesgue integral 


Stratonovich integral wrt q 


Stratonovich integral 


yields the following Stratonovich form of the stochastic EPDiff equation 

dm + £uixi = 0 . 


( 1 . 11 ) 


( 1 . 12 ) 


The momentum density m{x,t) and velocity vector field u{x,t) in (1.12) are related by to = H and the stochastic 
vector field u(x,t) is given by 


l{x, t) = u{x, t) dt — ^ ^i{x) o dWi{t ). 


(1.13) 


Proof. The first step is to take the elementary variations of the action integral (1.9), to find 

Si 

Su : - - poq = 0, Sp : dq + £zq = 0, 5q : —dp + £~p = 0. 

Su 


(1.14) 


The first equation in (1.14) follows from the definition of the diamond operation (o) in (1.10). The second and 
third equations immediately follow from variations of the equivalent form of the action S{u^p,q) in equation (1.9) 
and integrations by parts with vanishing endpoint and boundary conditions. The governing equation for to will be 
recovered by using the result of the following Lemma. 


Lemma 4. Together, the three equations in (1.14) imply (1.12). 

Proof. For an arbitrary p G X, one computes the pairing 

{dm , p)^ = {dp o q + po dq , p)^ 


By equation ( |1.14D = (^{£~p) o q - p o P^q , p 
\P : ( £u£r) 4“ £r}£u^Q 


V 


= {p , a.duP qJ^ = -{poq, aduP ^ ^ 
= - (ads(pog), 77^ =-(£um,p 


(1.15) 


X 


X 
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where ad^?? = —[<£ 5 , £,,] = ad^.ry o dWi(t) is the stochastic adjoint action differential. Since rj G X 

was arbitrary, the last line completes the proof of the Lemma. In the last step we have used the fact that coadjoint 
action (ad*) is identical to Lie-derivative action {£) for vector fields acting on 1-form densities. □ 


The result of Lemmanow produces the m-equation in (1.12) of Theorem]^ This completes the proof. □ 
Remark 5 (Multiplicative noise in the 3D vector stochastic EPDiff equation). 

In 3D vector notation, the 1-form density m is expressed as, m = m ■ dx ® d^a; and equation (1.12 1 becomes 

dm (j)j{mu^) + mjVu^^dt + ^ (^dj[m^f{x)) + mj V^iix))odW,it)=0, (1.16) 


with m = 6ij6u. Importantly, the noise terms in (1.16) multiply both the solution and its gradient. The latter 
is not a common form for stochastic PDEs. In addition, both the spatial correlations ^i{x) and their derivatives 
\7^i{x) are involved. The effects of these noise terms on the singular solutions of stochastic EPDiff in one spatial 
dimension will be treated in Section [2] and its numerical solutions will be studied thereafter. 


1.3 Ito version of the stochastic EPDiff equation 


In the Ito version of the stochastic EPDiff equation, noise terms have zero mean, but additional drift terms arise. 
These drift terms are double Lie derivatives, which are diffusive, as shown in m for stochastic fluid dynamics. 


The corresponding ltd forms of the stochastic EPDiff equation in (1.12) and the second and third equations in 
(1.14) are found by using Ito’s formula to identify the quadratic covariation terms |34| as 




dq -(- £ud ^ ^ ’ 


(1.17) 




dt , 


where £uq = £uqdt — J2i £^iqdWi{t) is an Ito stochastic differential related to the stochastic Lie differential (1.8) 
(it should be noted that £u is not a Lie differential) and we have used [dWi{t),dWj{t)\ = Sijdt for Brownian motion 
to identify the quadratic covariation terms as drift terms. For more details about this sort of calculation in the 
geometric mechanics context, see m- 


1.4 Legendre transform to Stratonovich stochastic Lie-Poisson Hamilton equations 


Theorem 6 (Lie-Poisson representation of Stratonovich stochastic EP equations). 

The Stratonovich stochastic EP system in (1.12) may be written equivalently in terms of a standard semidirect 
product Lie-Poisson Hamiltonian structure m with a canonical Poisson bracket {q,p}, as 


dm 


—ad(. jTTi 

-{■)oq 

po{-) 


Sh/Sm 

( 

m 

dq 

= 

-£(.)q 

0 

1 


Sh/5q 


q 

. dp. 


1 

-1 

0 


6h/Sp 

[ 

.p. 



(1.18) 


where h{m,q,p) is a stochastic differential representing the Legendre transform of the stochastic Lagrangian and 
{ •) ■}lp denotes the Lie-Poisson bracket. 
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Proof. As usual, the Legendre transform of the stochastic Lagrangian determines the stochastic Hamiltonian and 
its variational derivatives. In a slight abuse of notation, we may write this Legendre transform as 

h{m, q,p) = {m , u) dt — £{u, q)dt — E (m, ir{x))y^odWi{t) 

" (1.19) 

=: H{m,q,p)dt - ^ (m , ° dWi{t ), 


where now we allow g-dependence in the Lagrangian £(u,q). Varying the stochastic Hamiltonian in (1.19) gives 

'Sh 

6p 


L \ \ 5h ^ 

dti(m,q,p} = [ dm, — j + (du, m - — ) dt + (— , dq ) + [ — , dp 


= (^6m , udt — ^ o dWi{f^ + (^u, m — dt + ^ , dq'^ dt + {0 , dp) . 

Consequently, the corresponding variational derivatives of the stochastic Hamiltonian are 


( 1 . 20 ) 


dh 

dm 


udt — (x) o dWi(t) = u, 

i 


dh 

dq 


— —dt and 
dq 



( 1 . 21 ) 


The resulting Lie-Poisson Hamiltonian form of the system 
is then given by 


of Stratonovich stochastic variational equations in (2.5) 


dm 


--':4c 

— ad(. 

-(.)og 

po(-) 


dh/dm 


-ad~TO - o g 

dq 

= 

-£(.)q 

0 

1 


dh/dq 

= 


. dp. 


1 

-1 

0 


dh/dp 


1 

1 

1_ 


Of course, the terms involving dh/dq vanish, when dijdq = 0, as for EPDiff. 


□ 


Remark 7. The matrix operator in (1.22) is the Hamiltonian operator for a standard semidirect product Lie- 
Poisson structure [TB] with the canonical Poisson bracket (two-cocycle) {q,p} between q and p. Similar Lie-Poisson 
structures also appear in the Hamiltonian formulation of the dynamics of complex fluids [Hi- 


Stochastic Lie-Poisson Hamiltonian dynamics has also been previously studied and developed in applications in 
spin dynamics in iE] and in image registration in uniiisiES]. 

Canonical stochastic Hamilton equations were introduced in Bismut [1] and were recently developed further in 
the context of geometric mechanics in Lazaro-Cami and Ortega 


2 Stochastic variational perturbations in one spatial dimension 

2.1 Singular peakon solutions of the EPDiff equations 

The EPDiff(il^) equation in the one-dimensional case when £{u) = ^HuH^i = ^ J u'^ + a'^u'^dx is called the 
Camassa-Holm (CH) equation for m = di/du = u — a^Uxx with positive constant namely [8], 

mt + {um)x + mUx = 0 with m = u — a^Uxx ■ (2-1) 

This equation has singular peakon solutions, given by 

de ^ ^ 

m{x,t) := — = u - a'^Uxx = '^Pa{t)d{x - qa{t)), so that u{x,t) ■.= '^pb{f)K{x - qb{t)), (2.2) 

a— 1 b—1 

where K{x — y) = exp(—jcc — y\/a) is the Green’s function for the Helmholtz operator 1 — The peaked shape 

of the velocity profile for each individual peakon solution of the CH equation u{x, t) := p{t) exp(—|a; — q{t)\/a) gives 
them their name. 
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Figure 1: Singular peakon solutions emerge from smooth initial conditions and form a finite dimensional solution 
set for the CH equation, EPDiff(i7^). The velocity prohle for each individual peakon has the peaked shape given 
by u{x, t) := p{t) exp(—|a; — q{t)\/a), which is the Green’s function for the Helmholtz operator. The main point to 
notice is that the distance between any two peaks never passes through zero. That is, the peakons keep their order, 
even after any number of overtaking collisions. (The taller peakons travel faster.) 


Peakons are emergent singular solutions which dominate the initial value problem, as shown in Figure An 
initially confined smooth velocity distribution will decompose into peakon solutions and, in fact, only peakon 
solutions. Substituting the (weak) solution Ansatz (2.2) into the CH equation (2.1) and integrating against a 
smooth test function yields the following dynamical equations for the 2N solution parameters qa{t) and Pa{t) 


dqg 

dt 


u{qa{t),t) and ^ 
at 


-Pa{t) 


du{qa{t),t) 

dqa 


(2.3) 


The system of equations for the peakon parameters comprises a completely integrable canonical Hamiltonian system, 
whose solutions determine the positions qa(t) and amplitudes Paif), for all N solitons, a = 1,..., A^, and also describe 
the dynamics of their multi body interactions, as shown, for example, in Figure 


As mentioned earlier, the objective of the remainder of the paper is to use the Stratonovich stochastic EPDiff 
Theorem to study the effects of introducing this type of stochasticity on the interactions of the peakon solutions of 
the CH equation with Stratonovich parametric stochastic deformation (P-SD) and canonical Hamiltonian stochastic 
deformation (CH-SD). The first step is to adapt Theorem to accommodate the peakon solutions. For this 
adaptation, the advection condition (1.4) used previously will be replaced by the definition of peakon velocity as 
the time derivative of peakon position, as in the first equation in (2.3). 


2.2 Singular momentum map version of the Stratonovich stochastic EPDiff equations 


Theorem 8 (Canonical Hamiltonian Stochastic Deformation (CH-SD) of EPDiff ). 
The action S{u,p,q) for the stochastic variational principle SS = 0 given by 


S{u,p,q) 


i{u) dt + ^ {pa , dqa - u{qa,t) dt) j - / ^hi{q,p) o dWi(t), 


leads to the following Stratonovich form of the stochastic EPDiff equation 

dm = —£um dt + E {m, hi{q,p)} o dWi(t), 

i 

dqa = u{qa,t) dt + ^ {qa , hi{q,p)} o dWi{t) , 
i 

dpa = -Pa{t) — {qa,t) dt+ [pa , K{q,p)] O dWi{t) , 

i 


(2.4) 


(2.5) 
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where the momentum density m and velocity u are given by 


se 


N 


N 


{x,t) :=— = '^PaS{x-qa{t)), and u{x,t) :='^pbK{x - qb{t)). 


( 2 . 6 ) 


a—1 


6=1 


Proof. As in the proof of Theorem the first step is to take the variations of the action integral (2.4), to find 


5i 


N 


- Qait)) = 0 , 


dh,. 


Sp: dqa-u{qa,t)dt-^^{q,p)odWi{t) =0, 

^dpa 

Sq : -dpa-Pa{t) — {qa,t)dt-^^{q,p)odWi{t) = 0, 


(2.7) 


UX ^Qa 

i 

after integrations by parts with vanishing endpoint and boundary conditions. The first variational equation captures 
the relation (2.6), and latter two equations in (2.7) produce the corresponding equations in (2.5). Substituting the 
latter two equations in (2.7) into the time derivative of the first one yields the first equation in ( |2.5| ). □ 

The particular choice of the functions hi{q,p) = Paf,i{qa) reproduces the results of Theorem]^ for param- 

eterised stochastic deformation (P-SD) of the peakon solutions. We summarise this observation in the following 
Corollary. 

Corollary 9. [P-SD is a special case of CH-SD for EPDiff] Given the set of diffusivities f,i{x), i = 1,..., M, let 
hi{q,p) = J2a=iPa.iii^a)- Then the momentum density m{x,t) satisfies the equation 


dm -I- Tjjm = 0, 

where the stochastic vector field u{x, t) is given by the P-SD formula, 

u{x, t) = u{x, t) dt + ^ £,i{x) o dWi{t). 


( 2 . 8 ) 


(2.9) 


Proof. Specialise to hi{q,p) = Pa'?*(9a) the first line of equation ( |2.5[ ) in Theorem]^ 


□ 


Remark 10 (Outlook: Comparing results for P-SD and CH-SD). In Sectionj^and Sectionj^we will investigate the 
effects of choosing between two slightly different types of stochastic potentials on the interaction of two peakons, 
N = 2, corresponding to P-SD and CH-SD. The two options are h[^\q,p) = ^^^=1 Pa^iila) and h[^\q,p) = 
PaPia{q), respectively, for i = 1,... ,M. These are both linear in the peakon momenta and in the simplest 
case they have constant coefficients. We will consider numerical simulations for two cases: h^i\q,p) = c(pi + P 2 ) 
(P-SD for M = 1) and h^i\q,p) = fdipi, h^\q,p) = P 2 P 2 (CH-SD for M = 2) with constants c,/3i,/32- Although 
these choices for and are very similar, they will produce quite different solution behaviour in our numerical 
simulations of peakon-peakon overtaking collisions in Section 


Remark 11 (Stratonovich stochastic EPDiff equations in one dimension). 


1. In one spatial dimension, equation ( |2.8[ ) becomes 

dm + {um^ 2mUx)dt m^ f,i{x) o dWi{t) -\- 2m f'iix) o dWi{f) = 0 . (2-10) 

i i 

Importantly, the multiplicative noise multiplies both the solution and its gradient. The latter is not a common 
form for stochastic PDEs. In addition, both the spatial correlations £,i{x) and their derivatives are 

involved. 
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2. The equations for dqa and dpa in (2.5) are stochastic canonical Hamiltonian equations (SCHEs) in the sense 
of Bismut HUS]. These equations for dqa and dpa may be rewritten as 


^ TT _ OL 


dpa 


dpa 


where the deterministic Hamiltonian is given by 


n U _ . 

dpa = -T^{q,p)dt-Y^ -^{q,p)odW,{t) 
oqa ^ 


( 2 . 11 ) 


dqa 


H{d,P) = \ '^PaPbK{qa - db) ■ 


( 2 . 12 ) 


a.b 


The stochastic canonical Hamilton equations (2.111 can also be obtained by extremising the phase-space action 
functional 

N pT M 

S[q{t),p{t)] = (^Paodqa-H{q,p)dt) - j '^h,{q,p) o dWi{t). (2.13) 

Jo ^ _1 ^ Jo 


This is the restriction of (2.4) to the submanifold defined by the Ansatz 


3. In the ltd version of stochastic canonical Hamiltonian equations, the noise terms have zero mean, but additional 
drift terms arise. These drift terms are double canonical Poisson brackets, which are diffusive |26j : 

dp: dqa = u{qa,t)dt-\-'^{qa,hi}dWi{t) + ^ hj}, hi}dt , 


ltd Noise for q 


ltd Drift for q 


du 1 

Sq : dpa =-Pa-^iqa,t)dt-tr'^{pa,hi}dWi{t) + - . 


ltd Noise for p 


ltd Drift for p 


The ltd stochastic dynamics of landmark points in the image registration problem discussed in Trouve and 
Vialard [33 is recovered when we choose hi{p, q) = —aqi, with i = 1, 2, 3, for g G and constant cr. In that 
particular case, the double bracket terms vanish. In the present study, we will take hi{p,q) as the two cases 
mentioned in Remark |10| and compare their effects on the dynamics of peakons with K{x—y) = eyy>{—\x—y\/a) 
and pulsons with K{x — y) = exp(—(x — yYjc?) in one spatial dimension. We have introduced the latter 
Gaussian shaped pulsons, in order to determine how sensitively the numerical results we shall discuss below 
depend on the jump in derivative of the velocity profile for peakons. 


3 The Fokker-Planck equation 


The stochastic process in (2.11) for {q(t),p{t)) can be described with the help of a transition density function 


p{t,q,p;q,p) which represents the probability density that the process, initially in the state {q,p), will reach the 


state {q,p) at time t. The transition density function satisfies the Fokker-Planck equation corresponding to (2.11) 
(see [ 13 ], US])' Let us examine the form of this equation in the case h^i\q,p) = Pipi, h^\q,p) = j52P2- In that 


case the noise in (2.11) is additive, and the Stratonovich and ltd calculus yield the same equations of motion. 


3.1 Single-pulson dynamics 

Consider a single pulson {N = 1) subject to one-dimensional (i.e., M = 1) Wiener process, with the stochastic 


potential h{q,p) = Pp, where /3 is a nonnegative real parameter. The stochastic Hamiltonian equations (2.11) take 
the form dq = pdt P o dW{t), dp = 0, which are easily solved by 


qpit) =q+pt-\- pw{t), 


Pp{t) =P, 


(3.1) 

















DD Holm & TM Tyranowski 


Stochastic Solitons 


26 January 2016 


11 


where {q,p) are the initial conditions. Note that the pulson/peakon retains its initial momentum/height p. We will 
use this solution as a reference for the convergence test in Section [5.6[ The corresponding Fokker-Planck equation 
takes the form 


^ ^ _ n 

dt '^^dq 2^ 


(3.2) 


with the initial condition p(0, q,p; q,p) = 6 {q — q)5(p — p). This advection-diffusion equation is easily solved with 
the help of the fundamental solution for the heat equation, and the solution yields 


Pl3(,t,q,p;q,p) = 


1 _ (g-g-pt)^ 

e 2/5^* S{p — p). 


/3y^2Tri 


(3.3) 


This solution means that the initial momentum p is preserved, which is consistent with (3.1). The position has a 


Gaussian distribution which widens with time, and whose maximum is advected with velocity p. 

3.2 Two-pulson dynamics 

The dynamics of two interacting pulsons has been thoroughly studied and possesses interesting features (see [12], 
ppjl. It is therefore intriguing to see how this dynamics is affected by the presence of noise. Consider N = 2 
pulsons subject to a two-dimensional (i.e., M = 2) Wiener process, with the stochastic potentials hi{q,p) = Pipi 
and h 2 {q,p) = /32P2, where q = ( 51 , 52 ), P = (pi,P 2 ), and /3i,/32 > 0. The corresponding Fokker-Planck equation 
takes the form 




d 

dp2 


[a4q,p)p] - 


t.d^P 

dql 


' ifS = 0 


dq, 


with the initial condition p(0, q,p] q,p) = 6 {qi — qi) 6 {pi — pi) + 6 {q 2 — q 2 )S{p 2 — P 2 ), where 


ai{q,p) =Pi +P 2 K{qi - 52 ), 

02 ( 5 , 13 ) =P2 +PiK{qi - 52 ), 


03 ( 5 , P) = -PiP 2 K'{qi - 52 ), 
04 ( 5 , P) = PiP2K'{qi - 52 ). 


(3.4) 


(3.5) 


Despite its relatively simple structure, it does not appear to be possible to solve this equation analytically. It is 
nevertheless an elementary exercise to verify that the function 


p(t,5i,52,pi,P2;5i,52,Pi,P2) = pi3i{t,qi,pi;qi,pi) + pp^{t,q2,P2-,q2,P2), 


(3.6) 


where pp^ is given by (3.3), satisfies (3.4) asymptotically as 51 — 52 —± 00 , assuming the Green’s function and its 
derivative decay in that limit. This simple observation gives us an intuition that stochastic pulsons should behave 
like individual particles when they are far from each other, just like in the deterministic case. In order to study the 
stochastic dynamics of the collision of pulsons, we need to resort to Monte Carlo simulations. 

In Sectionj^we discuss our numerical algorithm, and in Sectionj^we present the results of our numerical studies. 


3.3 Two-pulson dynamics with P-SD 

The stochastic interaction of two (or more) pulsons can be analysed explicitly when the stochastic potential has the 


form h{q,p) = j5{pi + P 2 ) for ,0 > 0 (P-SD for M = 1; see Remark 10). It is an elementary exercise to show that in 


that case the stochastic Hamiltonian equations (2.11) are solved by 


9a(i) = 5a {t) + PW{t), 


Pait) =Pa (t). 


0 = 1 , 2 , 


(3.7) 


where q^it) and p^(t) are the solutions of the deterministic system (|2.3|). We verify this numerically in Section 5.5 
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4 Stochastic variational integrator 


Given the variational structure of the problem we have formulated in Theorem]^ it is natural to employ variational 
integrators for numerical simulations. For an extensive review of variational integrators we refer the reader to Mars- 
den & West [28] and the references therein. Stochastic variational integrators were first introduced in Bou-Rabee 
& Owhadi [Sj. These integrators were derived for Lagrangian systems using the Hamilton-Pontryagin variational 
principle. In our case, however, we hnd it more convenient to stay on the Hamiltonian side and use the discrete 
variational Hamiltonian mechanics introduced in Fall & West [^. We combine the ideas of |5| and , and propose 


the following discretization of the phase-space action functional (2.13): 


K-l N 


= E E -EE 


K-l M 




AW, 


k 5 


(4.1) 


k—0 i—1 '' ^ k=0 m—1 

where At = T/K is the time step, denote the position and momentum at time tk = kAt, and AWJ^ ^ 

N{0, At) are independent normally distributed random variables for m = 1,..., M and fc = 0,... ,K—1. Let L{q, q) 
denote the Lagrangian related to H{q,p) via the standard Legendre transform q^ = dH/dpi. Then one can easily 


show that (4.1) is equivalent to the discretization 


K-l / N k 


k+l _ 


At 


k—0 m—1 ^ 


(4.2) 


k=0 '' 

of the Hamilton-Pontryagin principle used in l^. Omitting the details, ( |4.1| ) is obtained by computing the left 
discrete Hamiltonian H~{q^^^,p^) corresponding to the discrete Lagrangian Ld{q^,q^^0 = AtL{q^^^, {q^^^ — 
q^)/At). The interested reader is referred to |25| for more details. Extremizing ( |4.1[ ) with respect to and p* 
yields the following implicit stochastic variational integrator: 

M 


- qt 

At 


dH 


dhr, 




AW, 


k-l 


AWl 


m—1 

M 


dpi 


k f)TJ OL 

^ pk\ _ k-\-l 


2At 

AWH 


AW 


At 


dqi 


,P 


fe-i-i 


2At 


(4.3) 


for i = 1,..., N. Knowing {q^,p^) at time tk, the system above allows to solve for the position q^~^^ and momentum 
at the next time step. For increased computational efficiency, it is advisable to solve the first (nonlinear) 
equation for first, and then the second equation for 

Note that in (4.1) we used pf — qf) to approximate the Stratonovich integral Pi o dqi in (2.13), which 


means the numerical scheme (|4.3|) will not be convergent for general nonlinear stochastic potentials hi{q,p) (i.e.. 


multiplicative noise). If we used {pf Tp^"^ )(9i — qi)l2, the resulting integrator would be a two-step method, 

i.e., it would not be self-starting, and its geometric/symplectic properties would be in question. Nevertheless, for 


additive noise, i.e., when the stochastic potentials hi{q,p) are linear in their arguments, the integrator (4.3) is a 
simple modification of the Euler-Maruyama method and its convergence can be established using similar techniques 


(see 129]). The integrator ( ]4.3| ) has strong order of convergence 0.5, and weak order of convergence 1. We further 
verify this fact numerically in Section |5.6| 


The integrator (4.3) is symplectic, and preserves momentum maps corresponding to (discrete) symmetries of 


the discrete Hamiltonian—for instance, if H{q,p) and all hi{q,p) are translationally invariant, as in our simulations 
in Section]^ then the total momentum 'Yl!i=iPi is numerically preserved. The proof of these facts trivially follows 
from |5], keeping in mind that the momenta pt and velocities qi are related via the Legendre transform. 


5 Numerical experiments 

We performed numerical simulations of the rear-end collision of two pulsons for two different Green’s functions, 
namely K{qi — ^ 2 ) = and K{qi — 92 ) = In the latter case, the corresponding pulsons are 
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commonly called ‘peakons’. We investigated the initial conditions qi = 0, q 2 = 10, p 2 = 1 together with the 
following four initial values: pi = 8, pi = 4, pi = 2, pi = 1. That is, we varied the initial momentum of the 
faster pulson. We perturbed the slower pulson by introducing a one-dimensional Wiener process with the stochastic 
potential h{q, p) = j3p2 (this corresponds to = 0, /32 = /? in Section 3.2). The pulsons were initially well-separated, 
so their initial evolution was described by (3.6). The parameter (3 was varied in the range [0,6.5]. We used the 
time step At = 0.02, and for each choice of the parameters 50000 sample solutions were computed until the time 
T = 100. 


5.1 Sample paths and mean solutions 



t t t 


Figure 2: Example numerical sample paths for Gaussian pulsons for the simulations with = 4 and /3 = 4. The 
positions are depicted in the plots in the upper row, and the corresponding momenta are shown in the plots in the 
lower row. 



t t t 


Figure 3: Numerical mean paths for Gaussian pulsons for the simulations with pi = 4. Results for three example 
choices of the parameter /3 are presented: (3 = 1.5 {left), /3 = 2.5 {middle), and /3 = 4.5 {right). The positions are 
depicted in the plots in the upper row, and the corresponding momenta are shown in the plots in the lower row. 

Figure shows a few sample paths from the simulations of the interaction of Gaussian pulsons for the case with 
Pi = 4 and (3 = A. The simulations for pi = 8 and pi = 2, as well as the simulations for peakons, gave qualitatively 
similar results. The most striking feature is that the faster pulson/peakon may in fact cross the slower one. In 
the deterministic case one can show that the faster pulson can never pass the slower one—they just exchange their 
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t t t 


Figure 4: Example numerical sample paths for Gaussian pulsons for the simulations with pi = 1 and /3 = 5. The 
positions are depicted in the plots in the upper row, and the corresponding momenta are shown in the plots in the 
lower row. 



t t t 


Figure 5: Numerical mean paths for Gaussian pulsons for the simulations with pi = 1. Results for three example 
choices of the parameter fi are presented: fi = 0.5 {left), /? = 1 (middle), and /3 = 2 (right). The positions are 
depicted in the plots in the upper row, and the corresponding momenta are shown in the plots in the lower row. 


momenta. The proof relies on the fact that both the Hamiltonian and total momentum are preserved (see [12] . 
PD]!. In our case, however, the Hamiltonian (2.12) is not preserved due to the presence of the time-dependent noise 
(see Figure]^, which allows much richer dynamics of the interactions. This may find interesting applications in 
landmark matching—see the discussion in Section 


Looking at Figure we also note that our variational integrator exactly preserves the total momentum, as 
expected. Figure depicts the mean solution for Gaussian pulsons with the initial condition = 4 for different 
values of the noise intensity /3. We see that for small noise the mean solution resembles the deterministic one, but as 
the parameter j3 is increased, the mean solution represents two pulsons passing through each other with increasingly 
less interaction. We study the probability of crossing in more detail in Section [5^ 


We observed that pulsons may cross even when they have the same initial momentum (see Figure]^. In the 
deterministic case they would just propagate in the same direction, retaining their relative distance. Nevertheless, 
the mean solution (see Figure]^ does not show any crossing. 


5.2 Probability of crossing 

We studied in more detail the distance between the pulsons Aq(t) = q 2 (t) — qi(t) at the end of the simulation, that 
is, at time t = 100. Figure [^presents the experimental probability density function of Aq computed for Gaussian 
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t 



t 


Figure 6: Mean Hamiltonian for the simulations with pi 


4 for Gaussian pulsons {left) and peakons (right). 



Figure 7: Numerical probability density p of the distance Aq{t) = q 2 {t) — qi{t) at time t = 100 for Gaussian pulsons 
for the simulations with pi = 4 (cf. Figure]^. Results for three example choices of the parameter fl are presented: 
(3 = 1.5 (left), j3 = 2.5 (middle), and j3 = 4.5 (right). 



Figure 8: Numerical probability density p of the distance Aq(t) = q 2 (t) — qi(t) at time t = 100 for Gaussian pulsons 
for the simulations with pi = 1 (cf. Figure]^. Results for three example choices of the parameter /? are presented: 
/3 = 0.5 (left), ft = 1 (middle), and (3 = 2 (right). 
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Figure 9: The probability of crossing, that is, the probability that 92 ( 1 ) < qi{t) at time t = 100, as a function of 
the parameter /3 for Gaussian pulsons {top) and peakons (bottom). 

pulsons with pi = 4. The density appears to have two local maxima, with the global maximum shifting from 
positive to negative values of Ag as the noise intensity /3 is increased. The simulations for pi = 8 and pi = 2, as 
well as the simulations for peakons, gave qualitatively similar results. Figure depicts analogous results for the 
case of Gaussian pulsons with pi = 1. In this case the density function also has two local maxima, but the global 
maximum never shifts to negative values of Ag. The probability of crossing as a function of the noise intensity /3 
is depicted in Figure We see that this probability seems to approach unity for the simulations with pi > 1, and 
0.5 for pi = 1. 


5.3 First crossing time 



Figure 10: The blue histogram presents the experimental probability density of the first crossing time Tc for 
Gaussian pulsons for the simulations with = 4 and /3 = 2.5. More precisely, this is the conditional probability 
density given that Tc < 00 , i.e., assuming the pulsons do cross, the integral p{t) dr yields the probability that 
the first crossing occurs at time Tc G [a, 5]. The black line depicts the inverse Gaussian distribution (5.3) with the 
parameters (5.2). 


It might be of interest to investigate the earliest time when the pulsons cross Tc = inf{t > 0 : Ag(t) = 0}—let us 
call it the first crossing time (also known as the first exit time or the hitting time] see |13j . |23|1. Assume Tc = +00 
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/3 


Figure 11: The mean first crossing time E(Tc) as a function of the parameter fi for Gaussian pulsons {top) and 
peakons {bottom). More precisely, this is the conditional expectation E{Tc\Tc < oo) given that the pulsons do cross 
(i.e., Tc < oo). 


if no crossing occurs. We can approximate the probability density of Tc using the single-pulson solution (3.1) and 
the asymptotic property (3.6). The two pulsons are initially far from each other, so we have 


Aq{t) Ki Aq-{pi-p 2 )t +l3W{t), (5.1) 

that is, Aq{t) is approximated by a Brownian motion starting at Aq = Aq{0) > 0 with the drift pi —p 2 . Assuming 
that Pi > p 2 , the probability density of Tc is given by the inverse Gaussian distribution Tc ~ IG{pL,X) with the 
mean p, and shape parameter A, respectively. 


A* = 


Aq 


P1-P2' 


where the density function is 


A = 


Af 


(5.2) 


Pig{t) = 


27TT^ 


exp 


-X{t- p,)'^ 

2p'^T 


(5.3) 


An example of a conditional probability density function of Tc is depicted in Figure [TOl It shows very good agreement 
between the experimental and approximate theoretical distributions. The corresponding density functions for all 
other simulations look qualitatively similar, differing in the mean, variance, etc. When pi —> then p —> +00 


and the inverse Gaussian distribution (5.3) tends to the Levy distribution, whose mean is infinite. However, we 


observed that for the simulations with pi = P 2 = 1 the first crossing time is still well-approximated by the inverse 


Gaussian distribution with the shape parameter A as in (5.2) and mean p which seems to asymptotically decrease 


with the noise intensity 13 . The conditional mean crossing time E{Tc\Tc < 00) (given that a crossing occurs) as a 
function of the noise intensity /3 is depicted in Figure [m The mean crossing time shows minor variations for the 


simulations with pi > 1 and agrees well with (5.2), while it appears to asymptotically decrease for the simulations 
with Pi = 1. 
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Figure 12: Example numerical sample paths for Gaussian pulsons for the simulations with the initial conditions 
qi = 0, pi = 4, q 2 = 10, and p 2 = 1, and the stochastic potential (5.4) with the parameters /3 = 4 and 7 = 4. The 
positions are depicted in the plots in the upper row, and the corresponding momenta are shown in the plots in the 
lower row. 




Figure 13: Example numerical sample paths for Gaussian pulsons for the simulations with the initial conditions 
qi = 0, Pi = 4, q 2 = 10, and p 2 = 1, and the stochastic potential h{q,p) = j3{pi + P 2 ) with the parameter /3 = 4. 


5.4 Noise screening 


In the numerical experiments described above we observed that the presence of noise causes pulsons to cross with 
a non-zero probability. The functions qi{t), Pi(t), 92 (^) and P 2 {t) define a transformation of the real line through 


(2.6). In the deterministic case this transformation is a diffeomorphism, but not when noise is added, since the 


crossing of pulsons introduces topological changes in the image of the real line under this transformation. This 
may be of interest in image matching, as in |35| . when one would like to construct a deformation between two 
images which are not exactly diffeomorphic. However, with that application in mind, one may want to restrict the 
stochastic effects only to the situation when two pulsons get close to each other. This can be obtained by applying 
the stochastic potential 


= PP2e' 


. ( 92 - 91 )^ 


(5.4) 


The parameter /3 > 0 adjusts the noise intensity, just as before, while the parameter 7 > 0 controls the range over 
which the stochastic effects are non-negligible. We performed a few simulations with this stochastic potential. Since 


this potential is nonlinear, the integrator (|4.3|) is not applicable here. Instead, we used the stochastic symplectic 
midpoint rule (see 


]). A few sample paths are depicted in Figure 12 
translation-invariant, so the total momentum is preserved. 


Note that this stochastic potential is 
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Figure 14: Example numerical sample paths for Gaussian pulsons for the simulations with the initial conditions 
qi = 0, Pi = 4, 52 = 10, and p 2 = 1, and the stochastic potential h{q,p) = j3q2 with the parameter /3 = 4. The 
positions are depicted in the plots in the upper row, and the corresponding momenta are shown in the plots in the 
lower row. 




Figure 15: Dependence of the absolute and mean errors on the time step for the single Gaussian pulson (left), two 
point vortices {center) and Kubo oscillator {right). 


5.5 Restriction to parametric noise and additive noise in the momentum equation 

Interestingly, crossing of pulsons does not seem to occur for the case of parametric stochastic deformation with the 
restriction (pia{q) = ^i{qa) as in Gorollary]^ We ran numerical experiments for the potential h{q,p) = I3{pi +P 2 ), 
which has the form as in Corollary with ^{x) = /3, but observed no interpenetration (see Figure 13). This 
is consistent with our observation in Section |3.3| and the fact that pulsons never cross in the deterministic case. 
We also did not observe crossing when the stochastic potential is independent of p. For instance, we performed 
simulations with the potential h{q,p) = j3q2. Such a potential results in additive noise in the momentum equation 
in (2.111 only, as in [3S]. A few sample paths are depicted in Figure]^ Note that in this case the total momentum 
is not preserved, since h{q,p) is not translationally invariant. In many cases the pulsons would asymptotically 
approach each other, but never pass. We observed similar behavior for the (translationally invariant) potential 
Hq,p) = /3exp(-(5i - 52 ) V 7 ) with /?,7 > 0. 
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5.6 Convergence tests 

5.6.1 Single pulson 


In order to test the convergence of the numerical algorithm (4.3) we performed computations for N = 1 Gaussian 
pulson subject to one-dimensional (i.e., M = 1) Wiener process with the stochastic potential h{q,p) = Pp (cf. 
Section 3.1). Simulations with the initial conditions g = 0, p = 4 and the noise intensity /? = 4 were carried until 
the time T = 2 for a number of decreasing time steps At. In each case 50,000 sample paths were generated. Let 
ZAt{t) = {qAt{t),PAt{t)) denote the numerical solution. We used the exact solution ( |3.1| ) as a reference for computing 
the absolute error E{\zAt{T) — z^(T)|) and the mean error \E{zAt(T)) — E{zi3{T))\, where Zjs^t) = {qis{t),pp{t)). 
The dependence of these errors on the time step At is depicted in Figure We verified that our algorithm has 
strong order of convergence 0.5, and weak order of convergence 1. 


5.6.2 Two planar point vortices 

We performed a similar test for N = 2 planar point vortices subject to a one-dimensional (i.e., M = 1) Wiener 
process. The system is described by 


Hiq,p) = rir2log 


\ai (72^ 



h{q,p) = pXiqi + yaipi + P\ 2 q 2 + 70 ' 2 P 2 , (5.5) 


where ri,r 2 are the circulations of the vortices, Oi = i/jrTfsgnTi, \i = \/|r 7 f are scaling factors, / 3,7 are the noise 
intensities, and qi,Pi denote the x- and y-coordinate of the i-th vortex, respectively (see E], m)- Simulations 
for Ti = 2, r 2 = 1 with the initial conditions qi = aiRi , q 2 = o' 2 i? 2 ) Pi = p 2 = Oj where i?i = r 2 /(ri -|- r 2 ), 
i ?2 = —ri/(ri -I- r 2 ), and the noise intensities /3 = 7 = 0.5 were carried out until the time T = 6.4 for a number of 
decreasing time steps At. In each case 50000 sample paths were generated. We used the exact solution (see [TT] . 

m) 


qi{t) = (JiRi cosujt + yaiW{t), Pi{t) = XiRi sinojt — pXiW{t), 

q2{t) = a2R2 cosujt + ya2W{t), p2{t) = X2R2 sinojt - pX2W{t), 


(5.6) 

(5.7) 


where oj = (Ti -|-r2)/(27r), as a reference for computing the absolute and mean errors (see Figure [I^. We verified 
that our algorithm has strong order of convergence 0.5, and weak order of convergence 1. 


5.6.3 Kubo oscillator 


To demonstrate that the integrator (4.3) fails to converge for multiplicative noise, we performed computations for 
the Kubo oscillator, which is defined by H{q,p) = p^/2 -|- q‘^/2 and h{q,p) = /3(p^/2 -|- 9 ^/ 2 ), where /3 is the noise 
intensity (see [SO])- The exact solution is given by 


q{t) = psm{t -I- pW{t)) qcos{t + pW{t)), 


p{t) = pcos{t + pW{t)) — qsin{t + pW(t)). 


(5.8) 


A similar convergence test with g = 0, p = 4, /3 = 1, and T = 6.4 revealed that the integrator (4.3) failed to 
converge, although the errors remained bounded (see Figure [T^. 


6 Summary 

We have seen in Section that the finite-dimensional peakon solutions for the EPDiff partial differential equation 
in one spatial dimension persist under both parametric stochastic deformation (P-SD) and canonical Hamiltonian 
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stochastic deformations (CH-SD) of the EPDiff variational principle. Being both finite-dimensional and canonically 
Hamiltonian, the dynamics of the peakon solution set for EPDiff admits the entire range of CH-SD in the sense of 
IS US], which includes P-SD but can be more general. Therefore, the peakon solution set offers a finite-dimensional 
laboratory for comparing the effects of P-SD and CH-SD on the stochastically deformed EPDiff SPDE solution 
behaviour. In fact, as it turns out, the peakon solution set for EPDiff offers a particularly sensitive assessment of 
the effects of stochasticity on finite-dimensional solutions of SPDE. In Section]^ we took advantage of the flexibility 
of CH-SD to study stochastic peakon-peakon collisions in which noise was introduced into only one of the peakon 
position equations (rather than symmetrically into both of the canonical position equations, as occurs with P-SD), 
while at the same time not introducing any noise into either of the corresponding canonical momentum equations. 
The precision and flexibility of the CH-SD approach to stochastic peakon-peakon collision dynamics revealed that 
its asymmetric case with noise in only one canonical position equation allows the sollton-llke singular peakon and 
pulson solutions of EPDiff to interpenetrate and change order on the real line, although this is not possible for the 
diffeomorphlc flow represented by the solutions of the unperturbed deterministic EPDiff equation. This crossing of 
peakon paths was observed and its statistics were studied in detail for CH-SD In numerical experiments in Sectionj^ 
In contrast, crossing of peakon paths was not observed for the corresponding P-SD simulations in which the noise 
enters symmetrically In both position equations. Crossing of peakon paths was also not observed when stochasticity 
was added only In the canonical momentum equations, as studied in |35j . 

Thus, for the deterministic EPDiff, adding stochasticity of constant amplitude with either CH-SD of P-SD to 
a finite dimensional Invariant solution set has been found to produce different SDE solution behaviour. Here, the 
difference has introduced the possibility of a topological change in the order of points moving on a line in the CH-SD 
approach, while no such change in topology seems to be available via the P-SD subclass. One can also Imagine that 
changing the level of noise in the P-SD EPDiff SPDE could change the number of peaks or pulsons; a feature which 
would not have been available if the level of noise were changed after the reduction to a fixed A^-peakon solution 
sector. The intriguing idea of creation of singular EPDiff solutions by P-SD noise in the SPDE is under current 
investigation. 

The investigation of stochastic EPDiff in this paper has raised and illustrated a potentially important issue. 
The need for assessing the validity of approximating the stochastic solution behaviour of nonlinear SPDE by SDE 
obtained from adding noise to finite-dimensional projections (or discretisations) of the solutions is likely to be 
encountered quite often in many other circumstances and can be expected to be of frequent future concern. In 
particular, this issue is likely to occur in considerations of model error in stochastic data assimilation. For example, 
the endeavours of computational anatomy must face this Issue In the use of the singular solutions of EPDiff known 
as landmarks in the task of registration of noisy images [35] . This issue of the validity of stochastic deformations of 
finite-dimensional approximations of evolutionary PDE is a challenge for continuing research in P-SD of EPDiff, as 
well as In stochastic deformations of more general continuum equations, such as Euler’s equations for an Ideal fluid, or 
the Navier-Stokes equations for a viscous fluid. The present work has shown that the introduction of even constant 
stochasticity into the equations of motion for exact solutions (peakons and pulsons, or landmarks, for EPDiff) 
can produce unexpected changes In topology of the solution in one dimension. The corresponding introduction of 
stochasticity into the equations of motion for finite-dimensional approximations such as discretisation, or projections 
of the solutions of nonlinear evolutionary PDE may result in other surprises. 
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